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1  INTRODUCTION 

The  application  of  the  principles  of  fracture  mechanics  to  the  safe  and 
economic  design  of  engineering  structures,  which  may  develop  cracks  in  service, 
requires  that  the  stress  intensity  factors  of  the  cracks  be  known.  In  any 
analysis  a  complicated  structure  must  usually  be  modelled  as  a  three-dimensional 
body  in  which  several  boundaries  interact  with  the  crack  tip.  The  stress  inten¬ 
sity  factor  will  now  vary  with  position  on  the  crack  front  and  cannot  be  accu¬ 
rately  represented  by  a  two-dimensional  model.  It  is  therefore  necessary  to 
develop  efficient,  accurate  and  economic  methods  for  determining  stress  intensity 
factors  in  three-dimensional  cracked  bodies. 

Techniques  based  on  weight  functions  are  efficient  since  once  known  they 
can  be  used  to  find  the  stress  intensity  factor  for  any  arbitrary  loading  from 
a  simple  summation  or  integration.  Weight  functions  for  three-dimensional  sym¬ 
metrical  crack  problems  have  been  determined1  using  the  finite  element  method 
(FEM).  In  this  work  it  was  observed  that  for  three-dimensional  problems  the 
boundary  element  method  (BEM)  may  lead  to  faster  solutions  because  of  the 
reduction  in  the  number  of  elements  that  would  be  possible.  The  BEM  has  been 
used  to  obtain  stress  intensity  factors  for  a  semi-elliptical  crack  subjected 
to  a  polynomial  distribution  of  tractions  over  the  crack  surfaces. 

3 

Recently  It  has  been  shown  that  an  application  of  Betti* s  theorem  allows 
the  stress  intensity  factor,  for  arbitrary  applied  stresses,  to  be  derived 
efficiently  and  economically  from  certain  ’weight  functions*.  These  weight 
functions  are  derived  from  the  boundary  displacements  of  the  body,  resulting 
from  point  forces  applied  near  the  tip  of  the  crack;  the  crack-tip  region  is 

i + 

modelled  by  a  Bueckner  singular  field  .  This  approach  was  developed  for  two- 
dimensional  bodies  containing  cracks  where  it  was  shown  that  a  boundary  method 
has  considerable  advantages  over  domain  methods5;  these  advantages  will  be 
greater  in  modelling  three-dimensional  bodies.  The  use  of  a  Bueckner  singular 
field  in  the  determination  of  weight  functions  reduces  the  number  of  boundary 
integrations  needed  In  the  BEM  calculations. 

In  the  present  application  to  three-dimensional  cracked  bodies  a  singular 
displacement  field,  analogous  to  the  two-dimensional  Bueckner  field,  is  derived 
from  the  analysis6  for  a  planar,  rectilinear  crack  subjected  to  equal  and  opposite 
localised  forces  on  the  crack  surface.  The  singular  displacement  field  is  then 
prescribed  over  a  small  cylindrical  cavity  along  the  crack  front  and  the  BEM  is 
used  to  obtain  weight  functions.  The  advantages  of  the  method  are  demonstrated 
for  an  edge-cracked  rectangular  bar  by  calculating  the  stress  intensity  factors 


A 


are  calculated,  from  the  weight  functions,  at  points  on  the  crack  front  for 
uniform  normal  tractions  on  the  ends  of  the  bar.  Even  for  relatively  coarse 
meshes  the  results  are  in  good  agreement  with  those  determined  from  other 
methods • 

2  FORMULATION  OF  METHOD 

For  an  elastic  body  Betti’s  reciprocal  theorem  can  be  written  as 


Ju?t*ds  =  jtfu'ds  »  (1) 

s  s 

where  l* ,  v?  are  the  self-equilibriated  tractions  and  displacement  field, 

A 

respectively,  corresponding  to  one  state  (state  1)  of  loading  and  _t  ,  u.  are 
those  corresponding  to  another  independent  state  (state  2)  of  loading.  The 
choice  of  ( t_  )  and  (£  ,li  )  are  arbitrary  except  that  they  must  separately 
satisfy  the  equations  of  equilibrium,  the  compatibility  conditions  within  the 
body  and  the  boundary  conditions  on  the  body. 

It  is  assumed  that  state  1  corresponds  to  the  desired  solution  of  a  crack 
in  the  body  whose  surface  is  subjected  to  a  set  of  self-equilibriated  tractions 
_t*  .  The  leading  terms  in  the  general  expansion  of  displacements  near  the  crack 
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tip  can  be  obtained  from  Sih  and  Liebowitz  in  the  form 


uf  mlLi —(F)  iKi  cos(!)[(1  '  2v)  +  sin2(l)]+  Ku 


(1  +  v)  /  2r ' 
E 


(irt  |Ki  sin(!)[: 


2(1  -  v )  -  cos 


!(!)]- 


t  m  _  (I+v)  /  2r 


u„  =  2 


(F)  Km  sin(l)  * 


sin(|)[2(l  -  v)  +  cos2(|)]j  , 

cos(l)[(l  -  2v>  "  sin2(l)]j- 


£  e 

where  u^  (j  »  1,2,3)  are  the  respective  local  displacements  in  the  local  x* 

coordinate  directions  parallel  to  the  normal,  the  binormal  and  the  tangent  to 

the  crack.  That  is  the  displacement  vector  u  for  state  1  can  be  written 


£-1  £-£ 

ulxl  +  U2X2  +  U3X3  * 


(3) 
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l 

where  x.  are  the  local  unit  vectors  in  the  x,  directions  (j  =  1,2,3).  The 
3  3  £  % 

coordinates  (r,0)  are  local  cylindrical  polar  coordinates  in  the  x^  ,  x^  plane, 

l 

centred  on  the  crack  front  at  =  0  .  The  stress  intensity  factor 

(N  *  I, II, III)  corresponds  to  the  opening,  sliding  and  tearing  modes  respec¬ 
tively;  E  and  v  are  the  elastic  modulus  and  Poisson's  ratio  respectively. 
Since  the  only  loading  to  be  considered  in  state  2  is  on  the  crack  surface,  it 
is  convenient  to  simplify  equation  (2)  to  obtain  the  displacement  components  on 
the  crack  surfaces  9  =  ±ir  at  a  distance  x^  behind  the  crack  front.  These 
are  given  from  equation  (2)  by 


2(1  - 

E 


■ 


2(1  - 

E 


2(1  + 


-  • 


The  state  2  is  taken  to  be  that  corresponding  to  a  solution  of  a  body 

containing  a  crack  which  is  subjected  to  pairs  of  concentrated  symmetrical 

forces  located  on  the  crack  surfaces  near  the  crack  tip  as  shown  in  Fig  1.  The 

tractions  t*  are  represented  in  terms  of  point  forces  in  each  of  the  three 
~ 

orthogonal  directions  and  are  given  by  the  vector  ?_  *  P^x^  +  ^2X2  +  ^3X3  w^ere 
P.(j  *  1,2,3)  are  the  magnitudes  of  the  point  force  in  the  x,  directions. 


For  point  forces  located  at  x* 


0  the  traction  vector  t 


can  be  written  as 


/  -  «(*?  * 

where  6 (x^  +  cj  represents  the  Dirac  delta  function.  Substituting  equations  (3) 
and  (7)  in  equation  (1)  gives 


P3KIII 

P1KII  +  P2K!  +TT^V) 


2/7(1  -  v  yc 


J  ~  ~ 
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Thus  the  stress  intensity  factor  for  any  mode  (N  =  I, II, III)  can  be 

determined  independently  by  making  P.  corresponding  to  that  mode  non-zero,  and 

J  * 

all  other  zero.  The  displacements  u  will  be  those,  on  the  surface  S  , 

corresponding  to  P.  acting  on  the  crack  near  the  tip.  Following  Bueckner's 

4  ^  & 

suggestion  for  two-dimensional  fields,  the  components  of  u  are  called  weight 
functions,  since  the  integral  in  equation  (8)  represents  a  weighted  average  of 
the  tractions  over  the  surface  S  .  Because  u*  depends  only  on  the  displace¬ 
ments  on  the  surface  of  the  cracked  body,  it  follows  that  a  boundary  method  of 
analysis,  such  as  the  BEM,  should  be  a  more  efficient  procedure  than  a  domain 
method.  The  boundary  integral  equation  to  be  solved  is 


Cij(-  )uj(-  }  * 


f TjjCx.x' )Uj(x)ds(x) 
S 


J  Uj j(x,x' )tj(x)ds(x) 
S 


(9) 


where 

tively  on  the  boundary  S 

tal  solution  in  an  infinite  elastic  (uncracked)  body  and  c 
constant . 


u^ ,  tj  are  general  displacements  and  tractions  (i,j  ■ 

U^j  >  and  )  are  the 

ij 


1,2,3)  respec- 
Kelvin's  fundamen¬ 
ts  a  known 


In  order  to  determine  u*  it  is  necessary  to  solve  the  boundary  integral 
equation  (9)  for  the  cracked  body  subjected  to  state  2  loading;  for  which 
equation  (9)  becomes 


c1j(x')u*(x')  + 


I  Ti:jOL»2L'  )u*(x)ds(x) 

s 


/ 


Uij(x,x' )t*(x)ds(x) 


(10) 


where  the  surface  integral  on  the  right-hand  side  of  equation  (10)  is  over  a 


small  near-tip  region  S 


0 


only,  since  t*  are  zero  everywhere  except  on  Sq 


Thus  the  computing  time  to  evaluate  the  right-hand  side  integral  is  less  since 
the  number  of  elements  on  Sq  is  much  less  than  the  total  number  of  elements 
on  S  . 


3  STRAIGHT  FRONTED  CRACK:  MODE  I 

Equations  (8)  and  (9)  are  now  applied  to  the  determination  of  weight 
functions  and  opening  mode  stress  intensity  factors  for  a  straight-fronted 
planar  crack.  For  this  case  P^  »  P^  *  0  in  equation  (8).  The  traction  t* 
now  becomes  the  applied  force  P^  .  In  order  to  avoid  numerical  difficulties 
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due  to  the  resulting  singular  field,  this  force  is  modelled  by  taking  to 

be  a  cylindrical  cavity  of  radius  r^  centred  on  the  crack  front  with  a  displace¬ 
ment  field  Vj  prescribed  over  its  surface. 


This  avoids  numerical  difficulties  due  to  the  singular  field  resulting 
from  and,  if  r^  is  small  enough,  adequately  represents  the  presence  of 

the  force  near  the  crack  tip.  The  displacement  field  is  taken  to  be  that 

for  a  semi-infinite  planar  crack  subjected  to  forces  perpendicular  to  the  crack 
surface  and  has  been  derived  from  the  general  solutions  given  by  Kassir  and 
Sih6 .  The  displacement  field  on  a  cylindrical  surface  r  *  Tq  near  the  crack 
line  may  be  written,  using  the  equation  (3)  in  the  Appendix,  with 
Uq  =  )  ,  pQ  =  p/rQ  ,  zQ  «  x3/rQ  and  Ri  =  lim  P2 *  as 


Equation  (10)  is  solved  for  u*  on  ^  an<*  on  Sq  ^or  a 

subjected  to  the  displacement  conditions  on  the  cylindrical  cavity  given  by 

equation  (11).  The  opening  mode  stress  intensity  factor  is  given  from 
equation  (8)  by 


H 


1  f  * 

- - - TTz  ±  •  £ 

1  -  v)(*r  )3/2  J 


d  s 


(12) 


where  U*  •  u*/u^  *s  normalised  weight  function  for  the  body  and  t*  are 
the  applied  tractions  for  which  1C  is  to  be  determined. 
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4  SOLUTION  STABILITY 

It  is  necessary  to  choose  an  optimum  value  for  the  radius  r^  which  must 
be  sufficiently  small  to  represent  adequately  the  concentrated  forces  near  the 
crack  front  but  not  so  small  that  the  cavity  cannot  be  modelled  without  using  an 

excessive  number  of  boundary  elements  on  its  surface.  The  stability  of  the 

solution,  as  r^  varies,  is  studied  using  an  edge  cracked  bar  of  rectangular 
cross-section,  thickness  t  ,  width  b  and  total  length  2h  .  The  bar  contains 
an  edge  crack  of  uniform  depth  a  across  the  thickness  t  . 

The  element  distribution  on  the  surface  of  the  bar  is  shown  in  Fig  2,  for 

a  =  0.5b,  t  =  1.5b,  h  =  1.5b  .  There  are  a  total  of  68  elements,  46  on  the  rec- 

tanglar  surfaces  and  22  on  the  cylindrical  cavity.  The  elements  on  the  cavity 
were  clustered  in  the  vicinity  of  the  origin  of  the  singular  field  in  order  to 
model  accurately  the  rapid  variation  of  the  field  in  that  region.  Because  of 
the  symmetry  of  the  bar  about  the  crack-line,  only  half  of  the  configuration 
needs  to  be  modelled. 

Weight  functions  were  determined  from  the  solution  of  equation  (10).  The 
stress  intensity  factor  at  the  mid-point  on  the  crack  front  was  calculated  from 
equation  (12)  for  a  uniform  stress  applied  over  the  ends  of  the  bar  as  shown  in 
Fig  3.  The  variation  of  the  normalised  stress  intensity  factor  /  o/ira  with 
r^/t  for  0.003  <  r^/t  <  0.01  is  shown  in  Fig  4. 

Results  for  various  values  of  t/a  and  h/a  are  compared  with  the  known 

Q 

plane  strain  value  which  they  should  approach  for  large  t/a  .  It  can  be  seen 
that  the  present  results  differ  from  the  plane  strain  value  by  less  than  2%  for 
Tg/t  <  0.006,  and  for  any  given  configuration  the  maximum  variation  is  12%  over 
the  whole  range  of  r^/t  .  For  subsequent  calculations  a  value  of  r^/t  *  0.004 
was  chosen. 

5  VARIATION  OF  THE  STRESS  INTENSITY  FACTOR  ALONG  THE  CRACK  FRONT 

In  order  to  illustrate  the  versatility  of  the  method  developed  in  this 
paper,  it  is  applied  to  the  determination  of  stress  Intensity  factors  at  differ¬ 
ent  positions  along  the  crack  front  for  the  edge  cracked  bar  shown  in  Fig  3, 

with  the  same  dimensions  as  given  in  the  previous  section.  This  problem  has 

9  10 

also  been  solved  by  Raju  and  Newman  using  the  FEM  and  by  Mendelson  and  Alara 

using  the  method  of  lines  (MOL).  The  boundary  element  equation  (10)  is  solved 
for  the  weight  functions  with  the  coordinate  origin  of  the  singular  field  at 
x^/t  »  0*0,  0.15,  0.3,  0.4  and  0.45  •  These  weight  functions  are  used  in 
equation  (12)  with  t?  -  ox^  where  o  is  the  uniform  tensile  stress  acting 
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on  the  ends  of  the  bar.  The  resulting  normalised  stress  intensity  factor 
Kj  /  o/ira  is  shown  in  Fig  5  as  a  function  of  position  along  the  crack  front. 
These  results  are  compared  with  those  determined  by  Raju  and  Newman  :  other 
results10  show  similar  trends.  Also  shown  in  Fig  5  is  the  plane  strain  solution 
for  the  case  t  ®  <*  .  This  is  close  to  the  three-dimensional  solution  near  the 
centre  of  the  body  but  differs  significantly  in  the  region  x^/t  >  0.4  ,  where 
the  traction-free  boundary  conditions  on  the  faces  of  the  bar  affect  the  stress 
intensity  factor. 

The  boundary  element  results  were  determined  using  68  surface  elements 

9 

compared  with  432  solid  elements  in  the  finite  element  analysis  .  Thus  the  BEM 
coupled  with  a  singular  field,  as  developed  herein,  represents  a  substantial 
saving  in  data  preparation  and  computer  storage.  Furthermore  since  the  method 
produces  weight  functions,  the  stress  intensity  factors  along  the  crack  front 
can  be  determined  for  any  symmetrical  loading  on  the  body  without  re-solving  the 
boundary  element  equations. 

6  CONCLUSIONS 

A  three-dimensional  singular  field  has  been  derived  and  used  to  model 
point  forces  acting  at  a  crack  tip  in  a  boundary  element  analysis.  The  tip  is 
replaced  by  a  cylindrical  boundary  over  which  the  field  acts;  an  optimum  value 
is  suggested  for  the  radius  of  the  cylinder. 

This  field  combined  with  the  boundary  element  method  can  be  used  to 
determine  weight  functions  on  all  surfaces  of  a  body  simultaneously,  and  hence 
the  stress  intensity  factors  for  arbitrary  loadings  can  be  calculated. 

This  combination  is  more  efficient  than  the  standard  boundary  element 
procedure  because  it  requires  less  boundary  integrations:  it  also  requires  much 
less  data  preparation  and  storage  than  domain  methods.  Even  with  relatively 
coarse  meshes,  accurate  values  of  the  stress  intensity  factor  were  obtained  for 
a  cracked  bar. 
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Appendix 

DERIVATION  OF  SINGULAR  FIELD 


From  Ref  6  the  displacement  field  for  a  semi-infinite  planar  crack  subjec¬ 
ted  to  opening  forces  ±P^  perpendicular  to  the  crack  are  given  (see  Fig  1  for 
Pj  =  0,  j  -  1,3  and  x^  =  x  ^ ,  j  =  1,2,3)  by 


00 


where 


-2(1  -  v)  f  +  x„ 


3f 
3  x0 


(i 


oo 


+  x0 


(A-l) 


2pTT  p 


-1 

tan 


The  constants  p  and  v  are  the  shear  modulus  and  Poisson's  ratio  respec¬ 
tively,  (r,6)  are  the  polar  coordinates  in  the  x^  ,  x^  plane,  p  is  defined 
by 


2  2  2 
(x1  +  c)  +  x2  +  x3 


and  the  coordinate  system  x 
the  Xj  directions  parallel 
crack  front. 


j 


to 


j  *  1,2,3  has  an  origin  on  the  crack  line  with 
the  normal,  the  binormal  and  the  tangent  to  the 


In  the  limit  as  c  ♦  0 


/2pp2 


/ r  sin 
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where 


and  now 


lim| 

c+0 


2  2  2  2 
P  =  XL  +  x2  +  x3 


Writing 


where 


equations  (A-l)  become 
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The  functions  A  and  B  are  then  given  by 

i  f2  [(xi +  *2)  -xJ  . 

A  -  -  j  I  - r - - -  dx.  -  t 

/  /  2  2  \  /  2  2  2  \ 

J  (x  +  x  (x  +  x  +  x  J 
00  \  1  2/  V  1  2  3' 


dx^  -  2x^ 


f2  [(*! +  4)  -  xi3 

/  2  2  2  \ 
lXl  +  x2  +  x3  ) 


B  =  -  2x„ 


f  ,(X1  +  X2  )  xi] 

(x2i  +  x2  +  x2) 


2  2  2 

Substitution  of  xi  +  x2  *  1  enables  the  integrals  to  be  evaluated  and 
the  final  expression  for  the  singular  field  becomes 


jj/2 


r3/2  sin  9  sln(|)  r  +  4r  cos  _ 

27  p3  U  “  J 


.1  .  EI  f  *2/2<1  -  U)ri  5l0(l)  ,  r3/2  ,l„28  r»  *r(l  -  CO.  6)ll  ' 

2  <7  2  ~T1  77\  l'  »  J[  7 

y  2  L  p  2/2  p  sin(^)  J 


(A-3) 


T  RT  2/2  x~r3  2  sin  0  sin(-y) 

v,  - - —  <  (1  -  2v)B - - - - - 

J  li/2  p4 


where 


—  \- 
y»  L2 


2x^/R 


/I  *„r3/2 


sin  +  T  cos 


PoJ  x. 


3r  (1  -  cos  0)  cosQ-^ 


-J7J  (co.  ♦„  +  ^  .in  *0)l 

(.In  *0  -  co.  , 
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where 


0 


L  =  In 

r(l  +  cos  0)  +  R  -  2/2r7  cos  $q  cos^)' 

r(l  +  cos  9)  +  R  +  2/ 2Rr  cos  cos(|j 

T 


tan-1 £ 


2/2  Rr 


sin  <\>0  cos 


r(l  +  cos  0)  - 


(l)i 

R 


R 


P 
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Fig  2  Representative  boundary  element  distribution 


